### sdef (Ratio Association method) analysis for healthy control (HC) vs ENL comparison

### Datasets - GSE129033 & GSE74481

library(sdef)
library(readxl)
library(dplyr)
library(ggplot2)
library(AnnotationDbi)
library(org.Hs.eg.db)
library(WriteXLS)

GSE129033_healthy_vs_ENL <- read.csv("GSE129033_healthy_vs_ENL.csv", header = TRUE)
GSE74481_healthy_vs_ENL <- read.csv("GSE74481_healthy_vs_ENL.csv", header = TRUE)


tmp.healthy_vs_ENL <- list(GSE129033_healthy_vs_ENL, GSE74481_healthy_vs_ENL)
names(tmp.healthy_vs_ENL) <- c("GSE129033_healthy_vs_ENL", "GSE74481_healthy_vs_ENL")

mylist <- list()

for(i in names(tmp.healthy_vs_ENL)) {
        df <- tmp.healthy_vs_ENL[[i]]
        df <- df[!duplicated(df$ENTREZID), ]
        df <- df[!is.na(df$ENTREZID), ]
        mylist[[i]] <- df
}

rm(i, df)

data <- Reduce(function(x, y)
      merge.data.frame(x, y, all = TRUE, by = "ENTREZID"), mylist)

data <- data[, c(1, grep("P_value.*", colnames(data)))]
row.names(data) <- as.character(data$ENTREZID)
colnames(data)[2:5] <- names(mylist)
data <- data[, 2:5]

keep <- (apply(data, 1, function(x)sum(abs(is.na(x)))) <= 1)
table(keep)
data <- data[keep, ]

# Replacing missing p-values with 0.5
data[is.na(data)] <- 0.5

# Remove ribosomal genes 
ribosome <- read.table("ribosomal_proteins.txt", header = TRUE, sep = "\t", stringsAsFactors = FALSE)
head(rib, 5)
data <- data[!row.names(data) %in% as.character(ribosome$gene_id), ]

setwd("...\Dataset_comparisons\healthy_vs_ENL\sdef")

#association test
Th_healthy_vs_ENL <- ratio(data = data, pvalue = TRUE, interval = 0.1)

Rh_healthy_vs_ENL <- baymod(iter = 500, output.ratio = Th_healthy_vs_ENL)

MC_healthy_vs_ENL <- Tmc(iter = 500, output.ratio = Th_healthy_vs_ENL)

feat.names <- row.names(data)

feat.lists.Bayesian <- extractFeatures.R(output.ratio = Th_healthy_vs_ENL, output.bay = Rh_healthy_vs_ENL, feat.names = feat.names, h = 0.1)

createTable(output.ratio = Th_healthy_vs_ENL, output.bay = Rh_healthy_vs_ENL)

annotation_healthy_vs_ENL <- select(org.Hs.eg.db, keys = as.character(feat.lists.Bayesian$User$`h= 0.1`$Names), columns = c("SYMBOL", "GENENAME"), keytype = "ENTREZID")

results_healthy_vs_ENL <- merge.data.frame(feat.lists.Bayesian$User$`h= 0.1`, annotation_healthy_vs_ENL, by.x = "Names", by.y = "ENTREZID") 

#Getting adjusted-Pvalues

data <- Reduce(function(x, y)
                merge.data.frame(x, y, all = TRUE, by = "ENTREZID"),
               mylist)

data <- data[, c(1, grep("^adjusted_P_value|^logFC.*", colnames(data)))]
colnames(data)[-1] <- paste(rep(c("logFC", "adjusted_P_value"), times = 4), rep(names(mylist), each = 2), sep = "_")
row.names(data) <- as.character(data$ENTREZID)
results_healthy_vs_ENL <- merge.data.frame(data[annotation_healthy_vs_ENL$ENTREZID, ], annotation_healthy_vs_ENL, by.x = "row.names", by.y = "ENTREZID") 
names(results_healthy_vs_ENL)

# Median log2FC
results_healthy_vs_ENL$Median_Log2FC <- apply(results_healthy_vs_ENL[, grep("^logFC", colnames(results_healthy_vs_ENL))], MARGIN = 1, FUN = median)

# DEGs with FDR < 0.1
results_healthy_vs_ENL$genes <- apply(results_healthy_vs_ENL[, grep("^adjusted_P_value", colnames(results_healthy_vs_ENL))], 1, function(x) sum(abs(x < 0.1)))

# Export output

WriteXLS(x = results_healthy_vs_ENL, ExcelFileName = "Results_healthy_vs_ENL.xls")
